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Abstract 

A hybrid scheme between large-scale electronic structure calculations is developed and applied 
to nanocrystalline silicon with more than 10 5 atoms. Dynamical fracture processes are simulated 
under external loads in the [001] direction. We shows that the fracture propagates anisotropically 
on the (001) plane and reconstructed surfaces appear with asymmetric dimers. Step structures are 
formed in larger systems, which is understood as the beginning of a crossover between nanoscale 
and macroscale samples. 
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Silicon is one of ideally brittle materials and its fracture behavior is studied intensively, 
because we can obtain essentially dislocation-free single crystals. A pioneering theory of 
brittle fracture was given, in 1920's, by Griffith^ within a continuum theory, which is the 
foundation of the present understanding of brittle fractures^. The fracture in single crystals 
should be investigated also in atomistic pictures, on the point of how and why the fracture 
path is formed and propagates in the crystalline geometry. This point includes surface 
reconstruction processes. Since fracture is a thermal non-equilibrium process, the atomic 
structure on a cleavage surface can be different from that on equilibrium clean surfaces. For 
example, the easiest cleavage plane in macroscale samples of silicon is the (111) plane, in 
which the surface structure is not the ground state (7 x 7) structure but a metastable 2x1 
structure 3 -^. 

This letter is devoted to the atomistic fracture behaviors in nanocrystalline silicon, espe- 
cially, its possible difference from macroscale samples. Such a difference can be expected, 
as explained below; Now a typical atomistic length in silicon is defined as do = \/vq « 3 
A, where vq gives the volume per atom. The essence of the Griffith theory^ is the energy 
competition between the energy gain of the strain relaxation and the loss of the surface 
formation energy. The former energy is a volume term proportional to (length) 3 , while 
the latter energy is a surface term proportional to (length) 2 . As analogous to the theory 
of nucleation^, the dimensional analysis gives the critical crack length for the spontaneous 
fracture propagation. The critical crack length cq is given as 1,2 

c g « (!) 

with the stress a, the Young modulus E(ia 10 2 GPa) and the surface energy per unit area 
7. The value of 7 was estimated to be in the order of U/m 2 ^, which can be reduced to 
the bond breaking energy (70^ « leV) in the atomistic picture. In a recent experiment 
with macroscale samples^, the stress is a ~ 10 1 MPa and Eq. (0) gives a macroscale length 
[cq ~ 1 mm). Since the length cq is not dependent on the sample size L, the fracture 
behavior can be expected to be different from the above picture, in case that the sample 
size L is smaller than the critical length cq (L < cq). In this letter, we will see such a 
situation in nanocrystalline silicon, in which the numbers of atomic layers for these lengths 
(« cc/do, L/do) are not macroscale numbers. 

For atomistic fracture simulations of silicon crystals, a recent work of classical modelings^ 
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FIG. 1: The computational time for bulk silicon as the function of the number of atoms, up to 
1,423,909 atoms; The CPU time is measured for one time step in the molecular dynamics (MD) 
simulation. A tight-binding Hamiltonian is solved using the perturbative order-N method and the 
exact diagonalization. We use a standard workstation with single CPU and 2 GB RAMs. 

was done with 10 5 atoms. A more recent work£, however, pointed out the limited applica- 
bility of classical modelings and the importance of electronic structure calculations. On the 
other hand, there are several ab initio calculations with 10 2 atoms&£. Due to the system size 
of simulations, these investigations are limited in situations, such as the preparation of the 
initial cleavage plane in which the reconstructed surface structure is assumed. Therefore, 
large-scale electronic structure calculations are essential. 

So far, we have developed several order-N methods for large-scale electronic structure 
calculations 10 . The order-N method is the general name of methods in which the compu- 
tational cost is proportional to the system size (N). We have developed the variational 
and perturbative order-N methods based on generalized Wannier state a 10 ' 11 . The Wannier 
states {4>i} are localized and the index % denotes its localization center. The equation for 
wave functions is common between the two methods and is given by the one-body density 
matrix (p = J2T C I n t ne computational algorithm, the perturbative method 

is simpler than the variational method. Figure Q demonstrates our large-scale calculation 
among 10 2 — 10 6 atoms. 

With the above two methods, we now construct a novel hybrid scheme, in which the 
variational method is used only for the wave functions whose centers locate near fracture 
regions. The regions contain, typically, 4 x 10 4 electrons. Some of such wave functions 
change their character dynamically from the bulk (sp 3 bonding) states to surface ones, as 
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discussed later. The other wave functions, in bulk regions, keep the character of the bulk 
bonding states and can be obtained by the perturbative method. The wave functions {&} 
calculated by the two methods are commonly used for constructing the one-body density 
matrix p. Any physical quantity is expressed by the density matrix^ and is well defined in 
the present hybrid scheme. 

The present work is based on a transferable tight-binding Hamiltonian with s and p 
orbitals^ 2 -. It is used for several crystalline phases and non- crystalline phases, such as liquicU 2 - 
and surfaces 1 ^. Since the fracture is the formation of surfaces in a bulk region, the theory 
should reproduce the atomic structures both in bulk and surface phases, which is satisfied 
in the present Hamiltonian. The essence of the quantum mechanical freedoms is the fact 
that the sp 3 -hybridized bonds are formed in the bulk region, while are not on surfaces. To 
analyze the hybridization freedom, a parameter fv) is defined, for a wave function (f)j, as 

/W = £ | (^Is) | 2 , (2) 

where \Is) is the s orbital at the 7-th atom. For example, / s = 1/4 in an ideal sp 3 hybridized 
state. 

In this letter, we focus on the Si(OOl) surface, a standard template of the modern silicon 
technology. A characteristic feature in the Si(OOl) surface is the formation of asymmetric 
dimer o 14 i 15 . The asymmetric dimer is connected by a ! er' bonding state. Another state is 
localized on the 'up' atom, the dimerized atom near the vacuum region. This localized state 
is called V state, because the direction of its p components is nearly perpendicular to the 
dimer bond. Here an energy quantity is defined as 



A4 cov) = (^H^i) - [f^e s + (1 - / S W HJ . (3) 

A negative value of Ae^' corresponds to the energy gain of a covalent bonding. The 'a' 
state has the gain of Ae\ cov ^ ~ — 2eV, which mainly contributes to the dimerization energy 
(about — 2eV) 1 ^. The V state has much smaller Ae\ cav \ which is comparable to the energy 
difference between the asymmetric and symmetric dimers (the order of O.leV)^. 

The simulation details are as follows; The hybrid order-N scheme is used for systems 
with 10 4 atoms or more. In systems smaller than the above size, the variational method is 
used in the whole region. The samples are isolated tetragonal clusters, whose geometries are 
labeled with the number of atomic layers, such as n 100 x n 010 x n 001 or n 110 x nn x n 001 . As 



the boundary condition, the Wannier states at the sample surfaces are terminated by fixed 
sp 3 bonding states and are not reconstructed. The time step of the molecular dynamics is 
3 fs. The total kinetic energy is controlled to be that with 300 K by the Nose thermostat 
method^. The numerical accuracy is checked among bulk, surface and fracture properties, 
such as the elastic constants, the dimer formations on the clean (001) surface, and the critical 
stress for fracture. The last quantity is checked, with smaller samples, in comparison with 
the standard diagonalization. The calculated fracture propagation velocity is always in the 
same order of, but less than, the Rayleigh surface wave velocity (4.5km/s), as expected from 
the continuum theory 2 -. 

For fracture propagations, external loads in the [001] direction are imposed. During 
the simulations, the external loads can be dynamically controlled by the atoms on the 
sample surfaces in the z direction. These atoms are fixed or under artificial constant- 
velocity motions in the z axis. The velocity, typically 10~ 2 km/s, is much smaller than 
that of observed fracture propagation velocities (km/s). As a seed of fractures, a short 
range repulsive potential is imposed on one particular pair of atoms, as a defect bond. For 
smaller samples, the simulations begin without initial deformations. The fracture always 
occurs with the external loads in the order of a « 1 GPa, which corresponds to the strain 
energy of ad\ ~ O.leV per atom. For larger samples, the simulations begin with initial 
static deformations in the above magnitude of external loads. The length cq in Eq. ((TJ is 
calculated as cq ~ lOOnm, which is longer than the present sample sizes (L < 20nm). 

In results, a two-stage reconstruction process is commonly observed as the elementary 
process during successive bond breakings; In Fig. Efa), we monitor the one-electron energy 
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FIG. 2: (a) Elementary reconstruction process with the one-electron energy Ej and the weight 
on s orbitals f®. (b) An asymmetric dimer on a resultant crack. The black rod and black ball 
correspond to the V and V states, respectively. 
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FIG. 3: Snapshots of a fracture process in the (001) plane. The sample size is reioo x ^oio x n ooi = 
33 x 33 x 33 (4501 atoms). The time interval between two successive snapshots is 0.3 ps, except 
that between (f) and (g) (about 1.3 ps). A set of connected black rod and black ball corresponds 
to an asymmetric dimer, as in Fig. E^b) . The left-down area has not yet fractured. 

Ei = ((fii\H\(j)i) and the hybridization freedom /W of a Wannier state \<pi). Before the bond 
breaking (t < ps), the wave function is a bonding state in the bulk region, deformed due 
to the external load. At t ~ ps, a bond breaking occurs and the wave function \(f>i) loses the 
bonding character with rapid increase of the bond length. Then (0ps<t<0.2ps), a twofold 
coordinated surface atom appears, since another bond is broken almost simultaneously. The 
wave function \<pj) forms a loan pair state that is stabilized by the increase of /W (0.6 0.8). 
The corresponding energy gain can be estimated to be — 0.2 x (e p — e s ) ~ — 1.3eV, which 
explains the energy gain in the figure (e^ = — 2.7eV — > — 3.8eV). In other words, the bond 
breaking process is caused by the local electronic instability, that is, the energy competition 
between the loss of the bonding (transfer) energy and the gain due to the increase of the 
weight on the s orbitals (/ s ). Finally, after the thermal motions with a finite time (t ~ 0.4ps), 
a pair of twofold coordinated atoms forms an asymmetric dimer with a a bonding state 
The corresponding covalent-bonding energy, defined in Eq. ©, is Aef ov) w -1.9 eV. This 
energy explains the gain in the figure (e^ = — 3.8eV — > — 4.8eV) and the energy loss (about 
1.3eV) due to the decrease of /W (0.8 — ■> 0.6). This asymmetric dimer is preserved until 
the end of the simulation, during a couple of pico seconds. Figure Efb) is an example of 
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observed asymmetric dimers. 




FIG. 4: (a) Ideal diamond structure with colored bond sites. (b)(c) Geometry of resultant cracks 
in the (001) plane. The broken bond sites are plotted as colored rods in the ideal (crystalline) 
geometry. Rods within one layer are painted in a same color, as in (a). The layer of black rods 
contains the defect bond at its central area. Atoms are plotted as dots. The sample sizes in (b) and 
(c) are niioXn 1 j xn ooi =49x50x49 (30025 atoms) and 97x100x49 (118850 atoms) , respectively. 
In (c), only a central area (niio x nn = 58 x 60) of the sample is shown. Note that the length of 
n no = 50 atomic layers is about 10 nm. 

Figure E2 shows the fracture process of a cubic sample with 4501 atoms. Each Wannier 
state is classified from its weight distribution into a bonding or atomic orbital, which is 
shown as a rod or a ball in the figures, respectively. A black rod or ball corresponds to one 
in the layer that contains the defect bond. One almost flat (001) surface is being created 
with many asymmetric dimers. The surface contains, however, many twofold coordinated 
atoms that have two back bonds (white rods) and a lone pair orbital (black ball). This 
is because a lone pair states are metastable, as discussed above. In Fig. |31 an anisotropic 
bond-breaking propagation is seen in the [110] and [110] directions, especially in the early 
snapshots. In the [110] direction, the successive bond breakings propagate along the nearest 
neighbor bond sites, which forms a zigzag path, as the black rods in Fig. HJa). A bond 
breaking process drastically weakens the nearest neighbor bonds, due to the local electronic 
instability, observed in Fig. Efa). Therefore, the successive bond breakings propagate easily 
in the [110] direction. In the [110] direction, on the other hand, the bond-breaking paths 
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are not connected, as the red rods in Fig. Ufa). In this direction, the bond breakings are 
propagated through the local strain relaxation, not by the local electronic instability. As 
results, the bond breaking propagation along the nearest neighbor bond sites (in the [110] 
direction of the present surface) is faster than that in the perpendicular direction (in the 
[110] direction), due to the difference of the successive bond breaking mechanisms. Note 
that a flat (001) surface is also obtained by a similar simulation without the initial defect 
bond, in which the fracture begins at the sample edges. 

Figures 0] (b) and (c) show larger samples with step formations^. In the two cases, all 
the conditions are the same, except the sample sizes. To see the step structures clearly, 
the broken bond sites are shown as rods in the ideal crystalline geometry. The defect bond 
is located in the center of the drawn area. The anisotropic fracture propagation in one 
(001) plane increases the anisotropic strain energy^. The anisotropy originates from the 
inequivalence between the [110] and [110] directions within one (001) layer. Since the above 
inequivalence does not appear within two successive layers, a step formation between them 
will release the anisotropic strain energy. In Fig. 0] (b), a step is formed between the layer 
of black rods and that of red rods. In the [110] direction, the bond-breaking propagation 
reaches the sample surfaces without step formations. In the [110] directions, the bond- 
breakings propagate slower and a step is formed in the central area at an early period of the 
crack propagation. After that, the fracture propagates among the two atomic layers. Since 
the two layers are symmetrically equivalent, the resultant step formation path, is almost a 
line in the [100] or [010] directions, as the boundary of the fractured areas between the two 
layers. 

In Fig Etc), the largest sample in the present letter, the above line structure does not 
reach the sample surfaces but is canceled with additional step formations in complicated 
paths. The sample size dependence of the step structures is understood by the beginning 
of the crossover between nanoscale and macroscale samples; If the sample contains so many 
atoms, the geometry of the resultant crack will be almost circular, as in Fig. 0] (c), so as 
to minimize the anisotropic strain energy^. If not, the strain energy is accumulated only 
within the confined bulk region due to the finite sample size. The resultant fracture behavior 
is directly related to the anisotropic atomic structure of the cleaved surface, as in Fig.|U(b). 

Since the above mechanism of step formations is two-dimensional, the present samples 
may be nanoscale 'thin' samples. In larger or thicker samples, an expected fracture behavior 
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is the bending of the fracture plane into the (111) plane, the easiest cleavage plane in 
macroscale samples, which is the crossover in the present context. In an enough large sample, 
the fracture mode with the easiest cleavage plane will grow with no regard for sample shapes 
and details of conditions. Note that the dynamical simulation with 10 5 atoms is the practical 
limitation within a single CPU workstation and the program code with parallel computations 
is now being developed for simulations with larger samples. 

This letter shows a possible difference in fracture behaviors among nanoscale and 
macroscale silicon crystals. Its origin is the size dependence of the energy competition 
between bulk and surface regions. The electronic structures between the two regions are 
essentially different and can be described by the present method with the well-defined total 
energy. This energy competition is also inherent in other phenomena, such as crystal growth 
and self organizations, which may be candidates for applications of the present method. 

This work is supported by a Grant-in- Aid for COE Research 'Spin-Charge- Photon' and 
a Grant-in-Aid from the Japan Ministry of Education, Science, Sports and Culture. This 
work is also supported by 'Research and Development for Applying Advanced Computational 
Science and Technology' of Japan Science and Technology Corporation. 



1 A. A. Griffith, Philos. Trans. R. Soc. London, Ser. A 221, 163 (1920). 

2 See textbooks of fracture, such as B. Lawn, Fracture of Brittle Solids, 2nd ed., Cambridge 
University Press (1993). 

3 K. C. Pandey, Phys. Rev. Lett. 47, 1913 (1981). 

4 F. Ancilotto, et al, Phys. Rev. Lett. 65, 3148 (1990). 

5 For example, L. D. Landau and E. M. Lifshitz, Statistical Physics, 3rd ed. Part I, Pergamon 
Press, Oxford (1980). 

6 J. C. H. Spence, Y. M. Huang, and O. Sankey, Acta metall. mater. 41, 2815 (1993). 

7 R. Perez and P. Gumbsch, Phys. Rev. Lett. 84, 5347 (2000). 

8 T. Cramer, A. Wanner, and P. Gumbsch, Phys. Rev. Lett. 85, 788 (2000). 

9 D. Holland and M. Marder, Phys. Rev. Lett. 80, 746 (1998). 

T.Hoshi and T.Fujiwara, J. Phys. Soc. Jpn, 69, No.12, 3773 (2000). 

1 T.Hoshi and T.Fujiwara, Surf. Sci. 493, 659 (2001). 



9 



I. Kwon, et al, Phys. Rev. B 49, 7242 (1994). 

For example, C.-C. Fu, M. Weissmann, A. Saul, Surf. Sci. 494,119 (2001). 
D. J. Chadi, Phys. Rev. Lett. 43, 43 (1979). 

A Ramstad, G. Brocks, and P. J. Kelly, Phys. Rev. B 51, 14504 (1995). 
S. Nose, Mol. Phys. 52, 255 (1984). 

The steps in the (001) surface are categorized in four types (D. J. Chadi, Phys. Rev. Lett. 59, 
1691 (1987)). The present surfaces, however, contain unreconstructed domains and are different 
from the above ones. 

The elastic property of silicon crystal shows only a small anisotropy within (001) plane; the 
[110] and [110] directions are equivalent and the values of the Young modulus are different by 
only about 30 % in the [100] and [110] directions. 



10 



